Bluefish wind analysis

1 Introduction

This wind analysis was inspired by Munch and Conover (2000), who found a correlation between wind speed and coastwide YOY index. The wind measurements used are described below:

The NOAA environmental buoy array provided environmental data from 1970 to 1992. Buoys collected hourly records of windspeed and wind direction at a height of 10 m, and surface temperature at a depth of 1.5 m. Data were examined from a total of 34 buoys in fixed locations throughout the region from Cape Fear to the Gulf of Maine and extending from nearshore to beyond the shelf break. Buoy locations are reported in Munch (1997).

Spring-spawned YOY were apparently characterized as fish between 13-27cm in length in the autumn bottom trawl, but “A complete set of length-frequency distributions and cohort designations is presented in Munch (1997)”.

And the observed correlation is described:

A strong relationship existed between offshore wind and recruitment. Using north (n) and east (e) windstress components from NOAA buoys near the edge of the shelf of the Mid-Atlantic Bight, the model ln(c.p.u.e.) = -16.5 n - 12.8 e + 0.59 (p<0.05) explained 36% of the variance in spring cohort abundance in the SNE region

I am trying to obtain a copy of Munch (1997), which is a master’s thesis (Munch, S. B. 1997. Recruitment Dynamics of Bluefish, Pomatonus saltatrix, on the continental shelf from Cape Fear to Cape Cod, 1973–1995. M.S. Thesis. State University of New York at Stony Brook).

2 Determine relevant strata

2.1 What polygons to use?

Entire east coast:

Inshore (ish) strata:

All strata south of ~ Cape Cod

3 Daily data

3.1 Downloaded daily mean dataset (uwind, vwind) from NOAA PSL

Located in data-raw folder.

3.2 Extract mean wind by day

3.2.1 v wind

3.2.1.1 Read in wind data

# use ecopull function
devtools::install_github("kimberly-bastille/ecopull")
vwind <- ecopull::nc_to_raster(nc = here::here("data-raw/vwnd.10m.2021.nc"),
                               varname = "vwnd",
                               show_images = TRUE) 

vwind <- raster::rotate(vwind)
raster::plot(vwind)

3.2.1.2 Mask to relevant strata

m_vwind <- raster::mask(x = vwind, 
                        mask = new_shape[new_shape$STRATA > 2000 & 
                                           new_shape$STRATA < 8000,2])
mean(m_vwind@data@values, na.rm = TRUE)

## [1] -0.9522257

m_vwind[[1]]

## class      : RasterLayer 
## dimensions : 30, 20, 600  (nrow, ncol, ncell)
## resolution : 1.031519, 0.6498195  (x, y)
## extent     : -80.45845, -59.82808, 30.21661, 49.71119  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +lat_0=40 +lon_0=-77 +x_0=0 +y_0=0 +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
## source     : memory
## names      : X2021.01.01 
## values     : 0.747654, 8.843357  (min, max)
## time       : 2021-01-01

raster::plot(m_vwind)

3.2.1.2.1 Plot over map

Seems like the existing strata are too small, try masking to more strata.

3.2.1.3 Mask to different strata

m_vwind2 <- raster::mask(x = vwind, 
                        mask = large_geom)
mean(m_vwind2@data@values, na.rm = TRUE)

## [1] -0.9731909

m_vwind2[[1]]

## class      : RasterLayer 
## dimensions : 30, 20, 600  (nrow, ncol, ncell)
## resolution : 1.031519, 0.6498195  (x, y)
## extent     : -80.45845, -59.82808, 30.21661, 49.71119  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +lat_0=40 +lon_0=-77 +x_0=0 +y_0=0 +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
## source     : memory
## names      : X2021.01.01 
## values     : 0.1705055, 8.967381  (min, max)
## time       : 2021-01-01

raster::plot(m_vwind2)

3.2.1.3.1 Plot over map

This seems like a better geometry.

3.2.2 u wind

3.2.2.1 Read in wind data

# use ecopull function
devtools::install_github("kimberly-bastille/ecopull")
uwind <- ecopull::nc_to_raster(nc = here::here("data-raw/uwnd.10m.2021.nc"),
                               varname = "uwnd",
                               show_images = TRUE) 

uwind <- raster::rotate(uwind)
raster::plot(uwind)

3.2.2.2 Mask to strata

m_uwind2 <- raster::mask(x = uwind, 
                        mask = large_geom)
mean(m_uwind2@data@values, na.rm = TRUE)

## [1] 1.321539

m_uwind2[[1]]

## class      : RasterLayer 
## dimensions : 30, 20, 600  (nrow, ncol, ncell)
## resolution : 1.031519, 0.6498195  (x, y)
## extent     : -80.45845, -59.82808, 30.21661, 49.71119  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +lat_0=40 +lon_0=-77 +x_0=0 +y_0=0 +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
## source     : memory
## names      : X2021.01.01 
## values     : 1.038357, 7.00613  (min, max)
## time       : 2021-01-01

raster::plot(m_uwind2)

3.2.2.2.1 Plot over map

3.2.3 Calculate weekly mean u and v wind

3.2.3.1 v wind

Positive v wind = wind is blowing to the north

##         names          mean_wind           min_wind         max_wind
## 1 X2021.01.01   5.94620376586914 -0.750392913818359 9.04062271118164
## 2 X2021.01.02   1.06632357895374  -3.64184379577637 9.03784370422363
## 3 X2021.01.03   1.17210961023966    -1.497154712677   3.040931224823
## 4 X2021.01.04 0.0323045295476913  -7.53194999694824 6.83328437805176
## 5 X2021.01.05  -4.79456804394722  -9.63622283935547 3.95362019538879
## 6 X2021.01.06   1.77418196439743  -10.3823204040527 5.81982803344727

3.2.3.2 u wind

Positive u wind = wind is blowing to the east

##         names         mean_wind          min_wind           max_wind
## 1 X2021.01.01  2.00335678100586 -5.33078384399414   8.30398178100586
## 2 X2021.01.02   10.223110652566  1.14752328395844   15.7158823013306
## 3 X2021.01.03  3.88333773295085 -4.47270393371582   8.53803825378418
## 4 X2021.01.04 0.293308817545573 -3.77553558349609   3.88755035400391
## 5 X2021.01.05 -4.23279141743978 -7.20564460754395   2.69376945495605
## 6 X2021.01.06 -8.83518941402435 -11.9605121612549 -0.417544364929199

3.3 Calculate across and along shore winds by week

Angles/directions are probably not reliable yet - have to ground truth

week	v mean	u mean	wind magnitude	wind blowing to	coarse angle	coast angle	adjusted angle	new angle	new wind blowing to	longshore wind	crossshore wind	check magnitude
1	2.025	-0.504	2.087	northwest	1.327	0.785	1.815	2.600	up_towards	1.788	-1.076	TRUE
2	-0.186	2.745	2.751	southeast	0.068	0.785	6.216	0.718	up_away	2.073	1.809	TRUE
3	1.335	-3.561	3.803	northwest	0.359	0.785	2.783	3.568	down_towards	-3.462	-1.574	TRUE
4	-6.751	-5.111	8.467	southwest	0.923	0.785	4.064	4.850	down_away	-1.160	8.388	FALSE
5	-5.425	2.651	6.038	southeast	1.116	0.785	5.167	5.952	down_away	-5.711	1.962	TRUE
6	-7.279	-5.137	8.909	southwest	0.956	0.785	4.098	4.883	down_away	-1.515	8.779	FALSE

3.3.1 Try plotting to see if things look correct

4 Monthly data

4.1 Extract mean wind by month

4.1.1 v wind

4.1.1.1 Read in wind data

# use ecopull function
vwind <- ecopull::nc_to_raster(nc = here::here("data-raw/vwnd.10m.mon.mean.nc"),
                               varname = "vwnd",
                               show_images = TRUE) 

vwind <- raster::rotate(vwind)
raster::plot(vwind)

4.1.1.2 Mask to relevant strata

m_vwind <- raster::mask(x = vwind, 
                        mask = large_geom)
mean(m_vwind@data@values, na.rm = TRUE)

## [1] -0.1615039

m_vwind[[1]]

## class      : RasterLayer 
## dimensions : 30, 20, 600  (nrow, ncol, ncell)
## resolution : 1.031519, 0.6498195  (x, y)
## extent     : -80.45845, -59.82808, 30.21661, 49.71119  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +lat_0=40 +lon_0=-77 +x_0=0 +y_0=0 +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
## source     : memory
## names      : X1979.01.01 
## values     : -0.188582, 2.251501  (min, max)
## time       : 1979-01-01

raster::plot(m_vwind)

4.1.1.2.1 Plot over map

4.1.2 u wind

4.1.2.1 Read in wind data

# use ecopull function
uwind <- ecopull::nc_to_raster(nc = here::here("data-raw/uwnd.10m.mon.mean.nc"),
                               varname = "uwnd",
                               show_images = TRUE) 

uwind <- raster::rotate(uwind)
raster::plot(uwind)

4.1.2.2 Mask to relevant strata

m_uwind <- raster::mask(x = uwind, 
                        mask = large_geom)
mean(m_uwind@data@values, na.rm = TRUE)

## [1] 2.585091

m_uwind[[1]]

## class      : RasterLayer 
## dimensions : 30, 20, 600  (nrow, ncol, ncell)
## resolution : 1.031519, 0.6498195  (x, y)
## extent     : -80.45845, -59.82808, 30.21661, 49.71119  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +lat_0=40 +lon_0=-77 +x_0=0 +y_0=0 +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
## source     : memory
## names      : X1979.01.01 
## values     : 0.6709892, 1.908766  (min, max)
## time       : 1979-01-01

raster::plot(m_uwind)

4.1.2.2.1 Plot over map

4.1.3 Calculate monthly mean u and v wind

4.1.3.1 v wind

Positive v wind = wind is blowing to the north

##         names          mean_wind          min_wind           max_wind
## 1 X1979.01.01   1.74945317900429 -0.93473881483078    2.8818564414978
## 2 X1979.02.01  -2.44980944072207 -3.76906156539917 -0.422381818294525
## 3 X1979.03.01   1.27876439795519 -1.22186589241028   2.07900977134705
## 4 X1979.04.01  -1.72685616961836  -3.5530412197113   1.93556571006775
## 5 X1979.05.01 -0.634130119963347 -1.85161530971527  0.681272804737091
## 6 X1979.06.01   2.66009804243843  0.27544116973877   3.98488116264343

4.1.3.1.1 Test correlation with recruitment index from the 2019 assessment model

4.1.3.2 u wind

Positive u wind = wind is blowing to the east

##         names          mean_wind             min_wind         max_wind
## 1 X1979.01.01   1.25624870378194 -0.00740670785307884 2.44020557403564
## 2 X1979.02.01   1.91080342319949    -2.98821616172791 5.48137092590332
## 3 X1979.03.01   3.74361495176951    0.520142674446106 6.84650373458862
## 4 X1979.04.01 -0.592399209913953    -2.45218873023987 1.38007032871246
## 5 X1979.05.01  0.284347456588876    -2.11945486068726 1.32381010055542
## 6 X1979.06.01   3.30452072178324    0.290787488222122 4.67998027801514

4.1.3.2.1 Test correlation with recruitment index from the 2019 assessment model

4.1.4 Rotate

4.1.4.1 Determine rotation angle

Angle from north: 0.7853982 (pi / 4)

4.1.4.2 Do the rotation

names	u mean	Year	DOY	week	month	v mean	wind magnitude	ns	ew	wind blowing to	coarse angle	coast angle	adjusted angle	new angle	new wind blowing to	longshore wind	crossshore wind	check magnitude
1979.01.01	1.256	1979	1979-01-01	1	1	1.749	2.153	north	east	northeast	0.948	0.785	0.948	1.733	up_towards	0.349	-2.125	TRUE
1979.02.01	1.911	1979	1979-02-01	5	2	-2.450	3.107	south	east	southeast	0.908	0.785	5.375	6.160	down_away	-3.084	0.381	TRUE
1979.03.01	3.744	1979	1979-03-01	9	3	1.279	3.956	north	east	northeast	0.329	0.785	0.329	1.115	up_away	1.743	3.552	TRUE
1979.04.01	-0.592	1979	1979-04-01	13	4	-1.727	1.826	south	west	southwest	1.241	0.785	4.382	5.168	down_away	-0.803	1.640	TRUE
1979.05.01	0.284	1979	1979-05-01	18	5	-0.634	0.695	south	east	southeast	1.150	0.785	5.134	5.919	down_away	-0.649	0.247	TRUE
1979.06.01	3.305	1979	1979-06-01	22	6	2.660	4.242	north	east	northeast	0.678	0.785	0.678	1.463	up_away	0.456	4.218	FALSE

Alongshore wind (positive = going up the coast)

Across shore wind (positive = going away from the coast)

4.1.4.2.1 Test correlation with recruitment index from the 2019 assessment model

4.1.4.2.2 Try GLMs

Poisson:

## 
## Call:
## glm(formula = Recruitment ~ 1, family = "poisson", data = dat)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2804.1  -1400.2   -413.9    759.6   6739.2  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 1.765e+01  2.525e-05  699034   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 128598930  on 33  degrees of freedom
## Residual deviance: 128598930  on 33  degrees of freedom
## AIC: 128599593
## 
## Number of Fisher Scoring iterations: 4

## 
## Call:
## glm(formula = eqn, family = "poisson", data = dat)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2781.3  -1250.3   -362.4    561.4   6077.9  
## 
## Coefficients:
##                         Estimate Std. Error   z value Pr(>|z|)    
## (Intercept)            1.765e+01  3.574e-05 494010.71   <2e-16 ***
## crossshore_wind_April  1.878e-02  1.401e-05   1340.83   <2e-16 ***
## crossshore_wind_May    2.667e-02  1.680e-05   1586.89   <2e-16 ***
## longshore_wind_April  -3.858e-02  1.464e-05  -2635.51   <2e-16 ***
## longshore_wind_May     1.614e-03  2.056e-05     78.52   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 128598930  on 33  degrees of freedom
## Residual deviance: 118338253  on 29  degrees of freedom
## AIC: 118338924
## 
## Number of Fisher Scoring iterations: 4

##      df       AIC
## mod1  5 118338924
## mod0  1 128599593

Negative binomial:

##  Family: nbinom2  ( log )
## Formula:          Recruitment ~ 1
## Data: dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1211.3   1214.4   -603.7   1207.3       32 
## 
## 
## Dispersion parameter for nbinom2 family (): 13.1 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 17.64736    0.04743   372.1   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##  Family: nbinom2  ( log )
## Formula:          
## Recruitment ~ crossshore_wind_April + crossshore_wind_May + longshore_wind_April +  
##     longshore_wind_May
## Data: dat
## 
##      AIC      BIC   logLik deviance df.resid 
##   1216.2   1225.3   -602.1   1204.2       28 
## 
## 
## Dispersion parameter for nbinom2 family (): 14.3 
## 
## Conditional model:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)           17.6560819  0.0648247  272.37   <2e-16 ***
## crossshore_wind_April  0.0205207  0.0266745    0.77    0.442    
## crossshore_wind_May    0.0264663  0.0294718    0.90    0.369    
## longshore_wind_April  -0.0410752  0.0277588   -1.48    0.139    
## longshore_wind_May     0.0005751  0.0365914    0.02    0.987    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##      df      AIC
## mod0  2 1211.333
## mod1  6 1216.186